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A set of interatomic pair potentials is developed for CdS and ZnS crystals. We show that a simple energy 
function, which has been used to describe the properties of CdSe [J. Chem. Phys. 116, 258 (2002)], can be 
parametrized to accurately describe the lattice and elastic constants, and phonon dispersion relations of bulk 
CdS and ZnS in the wurtzite and rocksalt crystal structures. The predicted coexistence pressure of the wurtzite 
and rocksalt structures, as well as the equation of state are in good agreement with experimental observations. 
These new pair potentials enable the study of a wide range of processes in bulk and nanocrystalline II-VI 
semiconductor materials. 



I. INTRODUCTION 

Many important processes in solid state materials, like 
the melting transition!^ structural transformations}^] 
or diffusion of impurities and defects^ require atomistic 
resolution in space and time for a comprehensive under- 
standing of the underlying mechanism. Despite major 
advances in electron microscopy, 6 - experiments can only 
provide a coarse-grained view of such processes. Molecu- 
lar dynamics computer simulations can in principle pro- 
vide the necessary microscopic perspective, but their pre- 
dictive power depends on the reliability and feasibility of 
available models. 

Methods based on first principles that retain a descrip- 
tion of the electronic state of the system potentially offer 
the highest accuracy. Because of their high computa- 
tional demand, however, they are not currently suited 
for in-depth studies of systems involving more than a 
few hundred atoms, or time scales longer than a few tens 
of picoseconds. Classical interaction potentials, parame- 
terized to reproduce emergent properties of the modeled 
material, offer a compromise between accuracy and com- 
putational speed. Potentials of different functional form 
and complexity have been developed for materials with 
widely disparate chemical and physical properties, rang- 
ing from water^ to golcP and biopolymersP Depending 
on the properties studied, agreement with experiment is 
usually good, in some cases rivaling or even besting that 
of affordable ah initio methods. 

Sparked by a comprehensive experimental study of the 
structural c hang es occurring in CdSe nanocrystals un- 
der pressureJ^E^HII] a simple pair potential has been de- 
veloped^ and successfully applied to reveal the mech- 
anisms of structural rear rangements in both bulk and 
nanocrystalline CdSepE^HH] However, simulation stud- 
ies of processes in many other semiconductor materials, 
or in multi-component systems like core/shell crystals 
or seeded nanorods, have been precluded by the lack of 
available models. Here we present a set of model po- 
tentials for crystalline CdS and ZnS, two semiconductors 
with potential use in various light harvesting and opto- 



catalytic devices W>- The potentials are designed to repro- 
duce the bulk lattice and elastic constants of the relevant 
crystal structures, as well as phonon dispersion relations. 
They are specifically constructed to be compatible with 
each other and with the existing model for CdSe™ and 
therefore also enable simulations of mixtures of the three 
compounds. 

The paper is organized as follows: In Section |TTJ we 
discuss the construction of the pair potentials and specify 
their parameters. In Section [III| we apply the models to 
calculate bulk enthalpies as a function of pressure and the 
equations of state for CdSe, CdS and ZnS and compare 
the predictions to experimental results. Discussion and 
conclusions are given in Section |TV| 

II. THE PAIR POTENTIAL 

We use the simple model developed for CdSe^ as a 
template to also describe CdS and ZnS. The two-body 
interatomic potential consists of a long range Coulomb 
part and a short range part which is represented by a 
Lennard- Jones (LJ) form: 
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where the indexes i and j refer to Cd, Zn, S, and Se 
atoms. To facilitate transferability and reduce the num- 
ber of parameters, we use standard combining rules for 
interactions of unlike atom types, namely — 

and Oij = \ (pi + Oj ) . 

The parameters g^, e^, and <Tj were obtained by fit- 
ting the lattice and elastic constants, and phonon dis- 
persion relations of bulk CdS and ZnS in three crystal 
structures: wurtzite, zinc-blende and rocksalt. As an 
additional constraint, the energy difference between the 
wurtzite and rocksalt struc ture at zero pressure was fit- 
ted to ab initio calculations J2S122I This was done to ensure 
that the wurtzite structure is the more stable structure 
at low pressures. The fitting calculations were performed 
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Table I. Potential parameters defining the interatomic inter- 
actions in ZnS, CdS and CdSe. 



at K, although the experimental data were obtained at 
finite temperatures. 

To obtain transferable potentials and ensure neutrality 
of the modeled materials, we fixed the magnitude of all 
ion charges to that of the original CdSe model, i. e., 
\qi\ = 1.18 e. We then proceeded in the following way. To 
obtain a model for CdS compatible with that of CdSe, we 
relaxed the LJ parameters for sulfur, keeping parameters 
for Cd fixed. In the second step, to arrive at a model 
for ZnS, we likewise relaxed the LJ parameters for Zn, 
keeping S fixed. Thus, a total of four parameters {uzm 
05, tzm an d eg) was determined with the fitting. 

We used a relaxed fitting procedure similar to the one 
discussed in Refs. l28l and 1291 This procedure is substan- 
tially more expensive computationally than the conven- 
tional fitting. However, it allows a much higher quality of 
fitting which is required to properly reproduce the struc- 
tural properties of CdSe. In the relaxed fitting procedure, 
the error was defined on the residual of the structural and 
dynamical properties of the optimized configurations of 
the different crystal phases rather than on the experimen- 
tal observed structures. Namely, the configurations and 
the lattice constants of each crystal phase were quenched 
using the conjugate gradient algorithm, and the afore- 
mentioned properties were calculated and compared to 
the experimental values for the quenched structures. In 
all the results reported here we have used Ewald sums to 
evaluate the electrostatic interactions,^ with a partition- 
ing parameter between the two spaces chosen to minimize 
the computational effort.^ The Lennard- Jones part of 
the potential was cut at half the box length (rs 10 A). 

The final parameter values resulting from the fit (in- 
cluding those for CdSe from Ref. 16) are summarized in 
Table [T] Consistent with their negative charge, the van 
der Waals (vdW) radii of the anions S and Se are sig- 
nificantly larger than those of the cations. In agreement 
with the corresponding ionic radii, the vdW radius of S 
is smaller than that of Se and the radius of Zn is smaller 
than that of Cd. For Zn we find that the best fit is 
obtained for a nearly vanishing vdW radius and a large 
value of e. These particular values are a result of the 
constraints imposed on the charge of the ions and the 
sequence of fitting steps. Because of the small value of 
a, the forces between Zn atoms at relevant distances are 
governed by the Coulomb term only; the large value of e 
sets the Zn-S bond length. A plot of the interatomic pair 
potentials is shown in Fig. [l] 

The accuracy of the model in reproducing the phonon 
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Figure 1. Interatomic pair potentials for Zn-Zn, S-S, and Zn- 
S (top panel) and Cd-Cd, Se-Se and Cd-Se (bottom panel). 
Open circles in the upper panel show the pure Coulomb re- 
pulsion term. 



dispersion relations of wurtzite ZnS, CdS and CdSe along 
the TA direction is shown in Fig. [2] The experimental 
results were obtained from inelastic neutron scattering^ 
on 116 CdSe and Raman scattering for ZnSp^l For CdS, 
we have compared our results with calculations provid- 
ing good agreement with infrared absorption.^ The split 
between the E2 and Bi; at the F-point indicates the ionic 
nature of the material.^ The overall frequencies in ZnS 
are higher due to the lighter masses of the atoms com- 
pared to CdS and CdSe. The TA and LA branches are 
back-folded into the lower E2 and Bi; branches, respec- 
tively. Similar back-folding occurs for the A l5 B^ and 
the upper E2 and Ei branches. As can be clearly seen in 
the figure, our simplified model captures the back-folding 
in all branches. The overall agreement between the cal- 
culations and the experimental results is reasonable given 
the simple form of the potential (cf. Eq. (JlJ). The model 
performs slightly better for ZnS and is more accurate for 
the lower frequency branches. The agreement can likely 
be improved usin g a p olarizable model, as is well known 
for alkali halidesP^SI We note that even better agree- 
ment wi th ex periment has been achieved with ab initio 
methods^ 21221 . 

In Table [TT] we compare the lattice and elastic con- 
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Figure 2. Phonon dispersion relations of wurtzite CdS (upper 
panel), ZnS (middle panel), and CdSe (lower panel) along 
the TA direction. The filled diamonds represent literature 
results™! 



stants calculated using our model with the correspond- 
ing experimental values. Note that the results for CdSe 
differ slightly from our original report.^ This deviation 
is a consequence of using a smaller tolerance for the min- 
imization, which was considered excessively cumbersome 
at the time the CdSe parameters were generated. The 
calculated lattice constants are within 1% of the experi- 
mental values, except for the case of the rocksalt struc- 
ture in ZnS, where the error is 3%. The agreement be- 
tween the calculated elastic constants and the experimen- 
tal values is qualitative, with small errors in Cn and C44. 

To compare the bulk modulus obtained from our model 
with experimental results we have used the well known 
relation between the bulk modulus and the elastic con- 
stants: 
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for hexagonal symmetry and 

B = [(C n + C22 + C 33 )/3 + 2C 12 ] /3, 
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for cubic symmetry. The calculated bulk moduli for 
wurtzite ZnS, CdS and CdSe are in reasonable agreement 
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Table II. Calculated lattice constants (in A) , elastic constants 
(in GPa), and bulk modulus (in GPa) of CdS, ZnS and CdSe 
for three crystal structures. Experimental results are given 
in parenthesis when available. (Elastic constants and bulk 
moduli were taken from Refs. l42l I44H461 1 



with the experimental values (which show some spread 
from one report to the other). Moreover, the values of 
the bulk modulus are compar able in accuracy to those 
obtained by ab-initio methods ! 26 l 27 l 



III. PHASE DIAGRAM AND EQUATION OF STATE 

To test the accuracy of our models on quantities not 
directly used in the fitting procedure, we calculated coex- 
istence pressures for the wurtzite and rocksalt structures 
at T = 300 K, as well as equations of state for all three 
crystal structures. 

In Fig. [3] we plot the enthalpy as a function of pressure 
for bulk ZnS, CdS, and CdSe. The pressure was varied 
between and 15 GPa using the constant pressure Monte 
Carlo simulation technique.^ For each crystal structure, 
we have used a periodically replicated simulation box 
of more than 400 atoms, and averaged the results over 
50,000 Monte Carlo cycles. Each cycle on average con- 
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Figure 3. Enthalpies as a function of pressure of bulk ZnS 
(upper panel), CdS (middle panel), and CdSe (lower panel), 
at a temperature of 300 K. The solid and dashed lines show 
results obtained for wurtzite and rocksalt crystal structures, 
respectively. The vertical dotted lines mark points of equal 
enthalpies and approximately correspond to thermodynamic 
wurtzite to rocksalt transition pressures in the respective ma- 
terials. 



sisted of one attempted displacement move for all atoms, 
and one attempted change of the simulation box volume. 

Approximating the true coexistence pressure by points 
of equal enthalpy, we find that the phase transforma- 
tion from wurtzite to the rocksalt structure occurs at rs 
2.4 GPa, w 1.6 GPa, and « 11.4 GPa for CdSe, CdS and 
ZnS, respectively. These values agree well with experi- 
mentally observed transition pressures of sa 2.5 GPap^EH 
« 2.5 - 3.2 GPapSffiD and « 12 GPa, 52 for CdSe, CdS, 
and ZnS, respectively. In all three materials, the zinc- 
blende crystal structure is not stable in the pressure range 
studied here; its enthalpy (not shown) is slightly higher 
than the corresponding wurtzite enthalpy. 

In Fig. [I] we plot the equation of state (volume as a 
function of pressure) for all three materials. We find 
excellent agreement with experiments on CdSe , 13 CdS,° 3 
and ZnSPl 



IV. CONCLUSIONS 

We have developed a set of transferable pair poten- 
tials for CdS and ZnS whose form is similar to that 
used for CdSeP^ The model consists of positively and 
negatively charged ions (Cd/Zn and S/Se, respectively) 
which interact via a Coulomb potential, supplemented 
by short-range repulsion terms and van der Waals at- 
tractive terms. In order to be able to model alloys and 
hetero-structures of CdSe/CdS/ZnS, we used standard 
combining rules for the cross terms and fixed the mag- 
nitude of the charges to the value obtained for CdSep^S 
thereby reducing the total number of fitting parameters 
to 4. The parameters were fitted to reproduce lattice 
and elastic constants, and phonon dispersion relations of 
wurtzite, zinc-blende and rocksalt crystals structures. 
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Figure 4. Volume as a function of pressure for bulk ZnS (up- 
per panel), CdS (middle panel), and CdSe (lower panel). The 
solid and dotted lines are the calculated results for wurtzite 
and rocksalt crystal structures, respectively. Vo is the unit 
cell volume of the wurtzite struc ture at zero pressure. Filled 
circles show experimental results.^ 3 * 52 ! 53 ! (Note that the CdSe 
data were obtained in experiments on 45 A diameter CdSe 
nanocr ystalsP) 



We have calculated the transition pressure of the 
wurtzite to rocksalt transformation, as well as equations 
of state at room temperature for all three materials. Our 
results are in good agreement with experiments, thus ver- 
ifying the accuracy and practicability of the pair poten- 
tial on quantities not used in the fitting procedure. 

As a final note, we point out that the simple functional 
form of the potential naturally limits the portability of 
our model. In particular, we have tested its accuracy 
in reproducing the lattice constants of ZnSe, a material 
whose properties were not included in the fit but which 
can be easily modeled using the parameters for Zn and 
Se. We find that deviations from experimental values are 
on the order of 5%, considerably larger than for the other 
materials. 

The current work extends our previous work on CdSe 
and provides a basis for the study of structural properties 
and dynamical processes of a larger variety of physically 
interesting materials. These include phase transforma- 
tion in core-shell structures, alloys which are important 
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for suppression of the Auger process ISilSSI seeded rods, 
and more. 
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